A fast impurity solver based on equations of motion and decoupling 
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In this paper a fast impurity solver is proposed for dynamical mean field theory (DMFT) based 
on a decoupling of the equations of motion for the impurity Greens function. The resulting integral 
equations are solved efficiently with a method based on genetic algorithms. The Hubbard and 
periodic Anderson models are studied with this impurity solver. The method describes the Mott 
metal insulator transition and works for a large range of parameters at finite temperature on the 
real frequency axis. This makes it useful for the exploration of real materials in the framework of 
LDA+DMFT. 
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I. INTRODUCTION 

Understanding exotic physical properties, such as high- 
T c superconductivity and the correlation-driven Mott 
metal-insulator transition of strongly correlated com- 
pounds (typically those including d or / electrons), re- 
mains a hard and fundamental task in modern condensed 
matter physics. During the past decade, the development 
and application of dynamical mean-field theory (DMFT) 
has led to a considerable improvement in our understand- 
ing of these systems The essence of DMFT is to map 
a many-electron system to a single impurity atom embed- 
ded in a self-consistently determined effective medium by 
neglecting all the spatial fluctuations of the self-energy 
However, this resulting quantum impurity model remains 
a fully interacting many-body problem that has to be 
solved, and the success of DMFT depends on the avail- 
ability of reliable methods for calculation of the local self- 
energy of the impurity model. 

Accordingly, much effort has been devoted to de- 
velop various impurity solvers. Among those, the iter- 
ated perturbation theory (IPT) 3 -! 4 -^, the non-crossing 
approximation (NCA) 7 -^^, equation-of-motion (EOM) 
methodiiii^ii^, Hubbard I approximation (HIA)- 15 , fluc- 
tuation exchange (FLEX) approximatio n 16 ! 17 , the quan- 
tum Monte Carlo method (Hirsch-Fye algorithm) (HF- 
QMC ) 18 ' 19 ' 20 ' 21 , the continuous time quantum Monte 
Carlo method (CTQMC) 22 ' 23 , the exact diagonalization 
(ED) 2 -i2£, the numerical renormalization group (NRG) 
metho d 26 ! 27 , and the density matrix renormalization 
group (DMRG) metho d 28 : 29 are widely adopted. How- 
ever, every impurity solver has its own limitation. IPT 
originally cannot be applied to the case away from half- 
filling while a modified IPT which can solve this prob- 
lem has to introduce an ansatz to interpolate the weak 
and strong coupling limits, and the generalization of IPT 
to the multi-orbital case requires more assumptions and 
approximations. NCA cannot yield the Fermi liquid be- 
havior at low energies and in the low temperature limit. 
The HIA can only be applied to strongly localized elec- 
tron systems like / electrons. FLEX works well in the 
metallic region while it fails in the large U region. Be- 



fore the appearance of CTQMC, the HF-QMC was not 
applicable in the low temperature limit and has serious 
difficulties in application to multi-orbital systems with 
spin-flip and pair-hopping terms of the exchange inter- 
action since the Hubbard-Stratonovich transformation 30 
cannot be performed in these systems. But even for 
CTQMC, the requirement to do analytical continuation 
of the results to the real frequency axis remains, which 
introduces some uncertainty especially for multi-orbital 
systems. In the ED method, an additional procedure 
is required for the discretization of the bath and as a 
consequence, the method is unable to resolve low-energy 
features at the Fermi level. NRG aims at a very pre- 
cise description of the low-frequency quasiparticle peaks 
associated with low-energy excitations while it has less 
precision in the Hubbard bands which are important in 
calculating the optical conductivity. Furthermore, all the 
numerically exact impurity solvers QMC, ED, NRG and 
DMRG are computationally expensive. 

However, today, a fast and reliable impurity solver 
is really urgently needed due to the fact that great 
achievements have been made in understanding cor- 
rectly the strongly correlated systems from first princi- 
ple by combining DMFT and local density approxima- 
tion (LDA) in density functional theory (DFT) , so called 
LDA+DMFT—. The aim of this paper is to present a fast 
and reliable impurity solver based on the EOM method. 
Equation of motion methods are limited by their decou- 
pling scheme, but EOM has shown its value by work- 
ing directly on the real frequency axis and at very low 
temperature. It can be a good candidate for a faster 
and reliable impurity solver by choosing a suitable de- 
coupling scheme. In fact, the infinite U case was studied 
by the EOM method for Hubbard model, periodic An- 
derson model and pd model in Ref. [ll|- In Ref. HH, the 
finite U case is studied without calculating the physical 
quantities self-consistently. Recently-i 3 -, this method has 
been improved by taking into account selfconsistency and 
applied to the Anderson impurity model in the large-iV 
limit. The operator projection method (OPM) 14 is re- 
lated to the EOM method. In this paper, we will use a 
different decoupling procedure than used previously for 
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a set of EOMs of the Anderson impurity model and then 
apply this new impurity solver to the finite U Hubbard 
model as well as periodic Anderson model via dynam- 
ical mean field theory. Meanwhile, we employ genetic 
algorithms to efficiently search for the self-consistent so- 
lution. The genetic algorithm significantly reduces the 
CPU time for convergence and improves the energy res- 
olution in the DMFT calculation. 

The paper is organized as follows: In Section [IT] we 
present the EOMs we use and introduce our decoupling- 
scheme. In Section IIIII we describe how the genetic al- 
gorithm is implemented in our DMFT loop. Finally, in 
Section HVl we test our EOM impurity solver on the Hub- 
bard model and periodic Anderson model. 

II. EQUATIONS OF MOTION AND 
DECOUPLING PROCEDURE 

We start with the Hamiltonian of the single impurity 
Anderson model. For arbitrary degeneracy N, it is given 
by 

+£( v & c L/*+wW (i) 

where c k(T , c kai fl and f a are the creation and annihi- 
lation operators for the conduction electrons and for the 
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where A(w) = ^2 k is the hybridization function 

and we have used 

(w-e k )((c k(T ;pJ = V k *AU-J i J. (10) 
These equations are generalized to arbitrary degeneracy 



correlated impurity electrons, respectively. n a = flf a 
corresponds to the density of the / electrons. e k is the 
dispersion of the conduction electrons, £/ is the site en- 
ergy of the correlated electron, U is the on-site Coulomb 
interaction strength of the / electrons, and V ka is hy- 
bridization between conduction and correlated electrons. 

In studying the system described by the Hamilto- 
nian of Eq. [T] we consider the double time temperature- 
dependent retarded Greens function in Zubarev nota- 
tion 32 , 

G AB (t, t') = ((A(t); B(t'))) = -<0(t - t')([A(t), B(t')} + )(2) 

involving the two Heisenberg operators A(t) and B(t'). It 
is convenient to work with the Fourier transform, which 
is defined as 

/oo 
dte^'-'VW; £(*')>>■ (3) 
-oo 

In the framework of the equation of motion method, the 
Greens function should satisfy the equation of motion 

u J ((A;B)) = ([A,B] + ) + (([A,H];B)) (4) 

where we have neglected the lower indices to. In the fol- 
lowing, all the Greens functions depend on frequency u>. 

As a result of the coupling between conduction and / 
electrons, we find the equations of motion, 
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N compared to Ref. [33|, i. e. they are at the same level 
as Ref. l35l Now a decoupling scheme is needed to trun- 
cate the equations of motion in order to get a closed set 
of equations. Here we have used the cluster expansion 
scheme proposed in Ref. where the higher order Greens 



i + (JV-iMW;4)U (5) 
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functions are separated into connected Greens functions 
of the same order and lower order Greens functions. The 
connected Greens function can not be decoupled any fur- 
ther as defined. This expansion scheme gives a natural 
and systematical way for truncation. It has been used in 
Ref. [33 for studying the single impurity Anderson model, 
in particular for infinite interaction strength U. This ap- 
proach to decoupling could be used to study the EOM 
method beyond the level of Ref. [H. The detailed cluster 
expansion scheme is given as 
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where digits 1-6 stand for operators, Asu t j >k -\ is the anti- 

symmetrization operator for operators (ijk), Spn^-,3,4) is 
the symmetrization operator for pair exchange between 
(1,2) and (3,4), and Greens functions or correlations 
marked by an index c represent connected terms. 

Using this decoupling scheme and neglecting all the 
three-particle connected Greens functions and those two- 
particle connected Greens functions which involve two c 
operators, i.e. ((c+ a ,c k , a , f a ; /]}} c , ((c^., f a . cy a \ /t)) c etc., 
and by assuming correlations with spin flip to be zero, 
e -9- (fl'fcr) = 0, we can get the single electron Greens 
function for arbitrary degeneracy N at the level of ap- 
proximation of Ref. [3J as 

(w - e d - A - AB){(f a ; /*» = 1 + A{n a , + C) (14) 



where 
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(17) 
(18) 



in which n a i — (ha'). This set of equations (fl4|) -(fT8l) is 
closed by the following two equations for the two-particle 
connected correlations: 



(vv) c = -J / d(uj')f(uj')lm((h a ,f a „-jl„)} c = -^ J ^O/^OM^/^;/^))-^^;/],,))) 
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where the two-particle Greens function can be obtained from the single-electron Greens function together with Eq. (J5]) 
and Eq. (J9j> . Finally the two connected correlations are 
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(21) 



Compared to Ref. [H, the set of equations (fbi ]) - (fl"5|) are 
generalized to arbitrary degeneracy AT. Ref. [35J has equa- 
tions of motion at the same level, but there the three par- 
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ticle Greens functions are neglected in the limit U — ► oo, 
while the Greens functions involving two c operators are 
considered to give little contribution for V — > 0. Thus, 
the decoupling of Ref. [H is constructed in the limit of 
parameters U — > oo, V — > 0. In Ref. [3?], the single im- 
purity Anderson model is studied for infinite interaction 
strength U and N = 2, and the focus is on the approxima- 
tion beyond that of Ref. [33J with the decoupling scheme of 
Ref. 3fl Here we have implemented the system of equa- 
tions (|14 p -(|18 p for finite U with arbitrary degeneracy N 
while neglecting the two particle connected correlations 
{h a „h a ,) c and (h a „cl al f a >) c . 

If we now specialize to degeneracy N — 2 and use 
the Hermitian conjugate (fl,Cka') = {c^'fa'), the set of 
equations (fT4"]) - (fT8|) becomes 
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with 
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We calculate the correlations self-consistently from the 
spectral theorem 

(flc ka ) = -- I duj'f(uj')lm V ^'f U;ftJ (25) 
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- y dw7(u/)Im{ 



<5fcfc' 
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where /(w') is the Fermi distribution function and the 
equation of motion for ((c| CT ; Ck'a)) has been used. 



method was employed to calculate the ground state wave 
function of one- and twodimensional quantum systems. 
We adopt this idea of optimizing the wave function un- 
til they obey the Schrodinger equation and carry it over 
to our optimization problem, that of finding a Greens 
function G(w) = that fulfils Eq. fl22J|. This ap- 

proach turns out to significantly improve the convergence 
speed, and as will be demonstrated below, the fact that 
it works with little or no broadening, the solutions are 
qualitatively better than from an iterative approach. The 
increase in convergence speed is essential for application 
of the SIAM model solution in DMFT calculations. 

The GA algorithm is started with a "population" of 
initial guesses. The imaginary parts of the initial pop- 
ulation of trial Greens functions are guessed as sums of 
Gaussians 



ImG(w) =L,(e 



+ e : 



2C^ 



(27) 



where L is a normalization factor, B, C are randomly gen- 
erated numbers and U is Coulomb interaction strength. 
We use the Kramers-Kronig relation to determine the 
real part 



ReGH = -i f l ^lcLo' 

7T ./ LO' — OJ 



(28) 



The convergence of the method can be speeded up if 
the positions of the randomly generated peaks cumulate 
around the positions of the Hubbard bands known from 
the atomic limit. Besides Gaussians, we have also tested 
other functional forms of the initial guess, but this had 
little influence on convergence speed and final converged 
result. 

The generation of trial Greens functions is evaluated 
and ordered according to a "fitness function" which mea- 
sures the closeness to a selfconsistent solution. Thus, we 
define the fitness function as 



III. METHODS OF SOLUTION 

In principle, the system of equations (f2"2"| - (|2"4"|) can be 
solved iteratively. But it turns out that the iterative 
solution requires significant Lorentzian broadening oj — > 
oj + ir/ and very small linear mixing factors a. Further- 
more, there are parameter regimes for which it is hard 
to converge a solution. The situation is not significantly 
improved by better mixing schemes like Broyden mix- 
ing^. Therefore, we turned to a different approach for 
finding the selfconsistent solutions. Genetic algorithms 
(GA) are adaptive heuristic search algorithms based on 
the idea of evolution by natural selectio n 39 i 40 and have 
been used in many optimization or minimization prob- 
lems of science and engineering^. In Refs.^HH, the GA 



F[G(u)] = \\G(u)-Tha[G(u)]\\ (29) 

where rhs[G(w)] represents the right hand sides of 
Eq. (fT4"l) or Eq. (|2"2"|) which are functional of G(oj) via 
the integral terms (JTEJ) - (HZJ) or ([25 ]) -([24 )> . The norm 

ll/HH = J <M/M| (so) 

measures the distance of the trial Greens functions from 
the selfconsistent solutions of Eqs. (fT4")l or l|22p . 

According to the standard procedure of GA, a new 
"generation" of the population of trial Greens functions 
is formed by application of two GA operators, "crossover" 
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Crossover: sample offspring 



Sample population members 




FIG. 1: (Color Online) Illustration for the genetic operations 
of crossover and mutation. The shown curves represent imag- 
inary parts of candidates for G{u). 



and "mutation" . The crossover operation is 
ImGf spring H = L{Im [G parent (u)O(lu - u a ) 

+ G p 2 arent (oj)e(uj -io)} 

ImGf spring (w) = L 2 lm \G parent (uj)®(ojq - u) 

(31) 

where loq is the randomly chosen crossover position, Q(lu) 
is the Heaviside function and L\, L 2 are normalization 
factors. Mutation introduces, with a low probability, ran- 
dom small changes in the trial Greens function in order 
to prevent the population from stabilizing in a local min- 
imum. The mutation operator is 



ImG offspring M = L(lmG parent (uj) + Ae 



(32) 



where A, B, C are randomly generated numbers and L 
normalizes the function. For both crossover and muta- 
tion, real parts are obtained via the Kramers-Kronig re- 
lation. Crossover and mutation operations are illustrated 
in Fig. [D 

Now the principles of selection have to be discussed. 
Some of the best members of a generation are preserved 
without replacing them by their offspring. Furthermore, 
the trial Greens functions that are actually included into 
the new generation are obtained by entering the result 
of the GA operations into Eqs. (fTl| or (|2"2"| and calcu- 
lating one iterative step. While in principle, pure GA 
operations could be used to find an optimal solution, the 
strict requirements imposed on a selfconsistent solution 
are more easily met by alternation of GA operations and 
iterative steps. As the new generation has more than 
twice as many members as the previous one, many are 
dropped according to their fitness values. In order to 
avoid premature convergence of the population to a sub- 
optimal solution, some members with unfavorable fitness 
values are kept in the population, and some new random 



trial Greens functions are added to the population. The 
end of the evolution is determined, as in the iterative 
solution of the integral equations, by a member of the 
population reaching the target accuracy. We usually use 
fitness function values of 10 -3 as a criterion for termi- 
nating the GA procedure. An additional advantage of 
the GA approach is the ease with which it can handle 
arbitrary kinds of constraints; they can be included as 
weighted components of the fitness function. 



IV. RESULTS AND DISCUSSIONS 

First, we investigated band- width control Mott metal- 
insulator transition in the Hubbard model. The densities 
of states (DOS) at four different values of U are shown 
in Fig. O As expected, quasi-particle peak as well as 
the upper and lower Hubbard bands are present in the 
metallic phase and transfer of spectral weight from quasi- 
particle peak to the Hubbard bands is clearly evident by 
reduction of the width of the central peak. In the insulat- 
ing state, the central peak suddenly vanishes and a gap 
appears between upper and lower Hubbard bands. Fur- 
ther increasing U leads to an increasing gap amplitude. 
The critical value of U for Mott transition obtained from 
our impurity solver is U c w 2.5. Compared to the criti- 
cal value from numerical renormalization group method 
where U c « 2.942£, our result underestimates the critical 
value of U due to the decoupling scheme. We note that 
in the metallic region, the height of our obtained DOS at 
the Fermi level is not fixed. This is due to the fact that 
two peaks in the imaginary part of the self-energy are 
quite close to the Fermi level, resulting in a numerical 
difficulty in getting a vanishing value of the imaginary 
part of the self-energy at Fermi level. 

Then, let us study the filling controlled Mott metal- 
insulator transition on the Hubbard model. In Fig. [31 we 
present the DOS as a function of doping at two different 
values of U. It is found that filling controlled metal- 
insulator transition occurs at 17 = 3 while the system 
remains in metallic state at U = 2. At U = 3, we also 
investigate the effective mass 
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<9Rc£(w) 



(33) 



as a function of doping concentration. It is shown in 
Fig. |4] that the effective mass clearly displays a diver- 
gent behavior as doping concentration goes to zero which 
seems to obey Brinkman-Rice picture for the Fermi liq- 
uid*^. In the small doping region, the carriers are more 
easily localized. We also studied the low frequency be- 
havior near the Fermi level for the metallic state at differ- 
ent temperatures, as shown in Fig. [5] We obtained that 
the imaginary part of the self-energy does not exactly 
follow Fermi liquid behavior under the present decou- 
pling scheme. However, as the temperature approaches 
zero, the negative imaginary part of the self-energy de- 
creases. The precision of the results at very low tern- 
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FIG. 2: (Color Online) DOS calculated with GA method for 
particle-hole symmetric Hubbard model on Bethe lattice 



perature is presently limited numerically. Therefore, the 
exact behavior of the imaginary part of the self-energy at 
the Fermi level at zero temperature is beyond our reach. 
Even though our decoupling scheme qualitatively shows 
an acceptable behavior, from principal considerations ex- 
act Fermi liquid behavior is not to be expected from a 
decoupling approach. 

We have also studied the Hubbard model with differ- 
ent types of bath as shown in Fig. [5] The influence of 
the bath has often been considered to be small since the 
self-consistent solution will not depend much on the ini- 
tial guess of the bath. Our result shows that both Bethe 
lattice and hypercubic lattice produce qualitatively sim- 
ilar results for the Mott transition. However, our results 
show that different baths yield different critical interac- 
tions strengths U c at which the Mott transition sets in. 
For the Bethe lattice, we find [/ c Bothc « 2.5, while for the 
hypercubic lattice, the result is £/ c hyporcubic « 2.4. 

We now turn to a comparison of the two methods of so- 
lution we employed, the iterative method with Lorentzian 
broadening and the combined genetic algorithm and iter- 
ative method. The comparison of the CPU time needed 
for a selfconsistent solution is clearly in favor of GA: 
While in general the solution with pure iteration takes 
four times longer, near the Mott transition the iterative 
solution becomes very slow and inefficient. In Fig. 
we show the results obtained with both methods for the 
same parameter values. In Fig. [JJ top left, we can see 
that in the DOS there is a nonzero continuous connection 
between two Hubbard bands in the result obtained with 
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FIG. 3: (Color Online) DOS for the asymmetrical Hubbard 
model on the Bethe lattice: filling controlled metal insulator 
transition. 




FIG. 4: (Color Online) Effective mass at different fillings for 
the Hubbard model on the Bethe lattice. 
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FIG. 5: (Color Online) Self-energy at low temperature close to 
the Fermi level for the metallic state (U — 1) for the particle- 
hole symmetric Hubbard model on the Bethe lattice. The 
inset shows the corresponding DOS. 
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FIG. 6: (Color Online) Comparison of Hubbard model with 
different kinds of bath: (left) Semicircle bath on Bethe lattice 
(right) Gaussian bath on hypercubic lattice, top two figures 
are metallic state, bottom two figures are insulator state just 
away from transition point 



Lorentzian broadening, which makes it difficult to distin- 
guish the Mott transition clearly when U approaches the 
critical value for the transition U c because the quasiparti- 
cle peak is very small in that case. Moreover, the Kondo 
peak will be greatly influenced by the amount of broaden- 
ing. Different broadening will give different critical value 
of U c . However, the combined GA and iteration method 
can give more precise results near the critical point which 
can be seen from the bottom left DOS figure. For the 
combined GA method, at U — 2.6 we find an insulat- 
ing state, while the Lorentzian broadening method still 
gives a metallic state at Fermi surface. This is due to the 
fact that the divergent behavior of the imaginary part of 
the self-energy just above the Mott transition cannot be 
correctly captured if there exists a finite broadening fac- 
tor. However, in the GA method, the broadening factor 
can be even set to zero, which eliminates the numerical 
problem induced by the factor. The right panels of Fig. [7] 
shows the comparison of the imaginary part of the self- 
energy. It is found that GA method really gives a correct 
divergent behavior even close to the Mott transition at 
the Fermi level, while the Lorentzian broadening method 
does less well. 

We have also studied the Hubbard model with arbi- 
trary degeneracy N. We have found that the decoupling 
scheme works nicely for N — 2, but for N > 2 there are 
some deviations from particle hole symmetry at half fill- 
ing. We observe that the band positions and occupation 
numbers are correct, but some broadening of the upper 
Hubbard band is missing. This shows that the presently 
used decoupling of the three-particle Greens functions 
misses some terms that would contribute to the damp- 
ing of the upper Hubbard band. We remedy this small 
deficiency by adding the ansatz S(N) = —c(N — 2)A(w) 
in the denominator of A (see Eq. IT51) because this de- 
nominator is mainly responsible for the upper Hubbard 
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FIG. 7: (Color Online) Comparison between the GA method 
and the iterative method with Lorentzian broadening for the 
particle-hole symmetric Hubbard model on the Bethe lattice: 
(top) DOS and self-energy for the insulating state with (7 = 3. 
(bottom) DOS and self-energy close to the Mott transition 
with U = 2.6. 



band, and neglected contributions from higher order 
Greens functions should contribute an unknown function 
of A(w). The factor c = 0.5 is found numerically from the 
requirement of particle hole symmetry at half filling, and 
the correction acts only for N > 2 as S(2) = 0. We get 
the results shown in Fig. where we have calculated the 
spectral functions for various degeneracies N at a tem- 
perature T = 0.01. With increasing N, the total on-site 
Coulomb interaction increases so that the two Hubbard 
bands shift further away from the Fermi level. But the 
critical U c also increases with N. Therefore, at the same 
U, the system shows more metallicity for larger N, and 
transfer of spectral weight is observed from upper and 
lower Hubbard bands to the Kondo peak with increas- 
ing N. This result is consistent with the QMC result of 
Ref. 51 

For comparison, we have calculated the Periodic An- 
derson model with our code in Fig. [9] We observe a 
similar behavior of the spectral weight transfer as in the 
large N Hubbard model. Meanwhile, Fig. [9] differs from 
the behavior shown in Ref. [ljl 



V. CONCLUSIONS 

We have presented the derivation and implementation 
of a solution of the single impurity Anderson model based 
on equations of motion and truncation. We employ a 
combination of genetic algorithms and iteration to solve 
the resulting integral equations. We demonstrate that 
our method is useful as an impurity solver in the context 
of dynamical mean field theory. We show results for the 
Mott metal insulator transition as a function of interac- 
tion strength U and as a function of filling n, and also 
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show the trend of weight transfer at different values of 
the spin-orbital degeneracy N. 



FIG. 8: DOS for the finite U Hubbard model for several values 
of the spin-orbital degeneracy N. The interaction strength is 
(7 = 2, the hybridization strength V — 0.25. (a) half-filling, 
(b) fixed impurity position Sf. 




FIG. 9: (Color Online) DOS for the periodic Anderson model 

for different values of the spin-orbital degeneracy N. (a) cor- We gratefully acknowledge support of the DFG 

related band, (b) conduction band. through the Emmy Noether program. 
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